ANALYSIS OF AN INTERFACE STABILISED FINITE ELEMENT 
METHOD: THE ADVECTION-DIFFUSION-REACTION EQUATION 

GARTH N. WELLS* 

Abstract. Analysis of an interface stabilised finite clement method for the scalar advection- 
diffusion-rcaction equation is presented. The method inherits attractive properties of both continuous 
and discontinuous Galerkin methods, namely the same number of global degrees of freedom as a con- 
tinuous Galcrkin method on a given mesh and the stability properties of discontinuous Galerkin 
methods for advection dominated problems. Simulations using the approach in other works demon- 
strated good stability properties with minimal numerical dissipation, and standard convergence rates 
for the lowest order elements were observed. In this work, stability of the formulation, in the form 
of an inf-sup condition for the hyperbolic limit and coercivity for the elliptic case, is proved, as is 
order k + 1/2 order convergence for the advection-dominated case and order k + 1 convergence for 
the diffusive limit in the L 2 norm. The analysis results are supported by a number of numerical 
experiments. 
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1. Introduction. Discontinuous Galerkin methods have proven effective and 
popular for classes of partial differential equations, in particular transport equa- 
tions in which advection is dominant. The attractive stability properties of suitably 
constructed discontinuous Galerkin methods and the possibility of matching non- 
conforming meshes are advantageous, but do come at the cost of an increased num- 
ber of global degrees of freedom on a given mesh compared to continuous Galerkin 
methods. In a number of recent works, advances have been made in reconciling the 
appealing features of continuous and discontinuous Galerkin methods in one frame- 
work. Works in this direction include those of Hughes et al. Labeur and Wells 2| 
and Cockburn et al. Q for the advection-diffusion equation, Burman and Stamm 4 1 
for advection-reaction equation, and Labeur and Wells Q and Labeur and Wells 5| 
for the incompressible Navier-Stokes equations. These methods generally strive for a 
reduction in the number of global degrees of freedom relative to a conventional dis- 
continuous Galerkin method without sacrificing other desirable features. In this work, 
stability and convergence estimates are presented for one such method applied to the 
scalar advection-diffusion-reaction equation, namely the interface stabilised method 
as formulated in Labeur and Wells [2|. 

The principle behind the interface stabilised method is simple: the equation of 
interest is posed cell-wise subject to weakly imposed Dirichlet boundary conditions 
in the spirit of discontinuous Galerkin methods. The boundary condition which is 
weakly satisfied is provided by an 'interface' function that lives only on cell facets 
and is single-valued on cell facets. An equation for this additional field is furnished 
by insisting upon weak continuity of the so-called 'numerical flux' across cell facets. 
This weak continuity of the numerical flux is in contrast with typical discontinuous 
Galerkin methods which satisfy continuity of the numerical flux across cell facets 
point-wise by construction. For particular choices in the method, it may be possible 
to achieve point-wise continuity. Upwinding of the advective flux at interfaces can 
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be incorporated naturally in the definition of the numerical flux, as is typical for 
discontinuous Galerkin methods. By building a degree of continuity into the interface 
function spaces (at cell vertices in two dimensions and across cell edges in three 
dimensions), the number of global degrees of freedom is equal to that for a continuous 
Galerkin method on the same mesh. The key to this reduction in the number of global 
degrees of freedom is that functions which are defined on cells are not linked directly 
across cell facets, rather they communicate only via the interface function. Therefore, 
functions on cells can be eliminated locally (cell-wise) in favour of the functions that 
live on cell facets. Outwardly the approach appears to have elements in common with 
mortar methods, and could serve to elucidate links between mortar and discontinuous 
Galerkin methods. 

The motivation for analysing the interface stabilised method comes from the ob- 
served performance of the method for the advection-diffusion in Hughes et al. [l| and 
Labeur and Wells 0] and for the incompressible Navier-Stokes equations in Labeur 
and Wells 0], and for the Navier-Stokes equations on moving domains, as presented 
in Labeur and Wells The method was observed in simulations to be robust and 
only minimal numerical dissipation could be detected. Labeur and Wells [2| also 
showed that the methodology can lead to a stable formulation for Stokes equation 
using equal-order Lagrange basis functions for the velocity and the pressure. The 
method examined in this work is closely related to that formulated by Hughes et al. 

1 1 for the advection-diffusion equation, and analysed in Buffa et al. [6| . Buffa et al. 

6| proved stability for a streamline-diffusion stabilised variant of the method, but not 
for the original formulation. For the case without the additional streamline diffusion 
term, stability was demonstrated for some computed examples by evaluating the inf- 
sup condition numerically. However, in the absence of an analytical stability estimate 
convergence estimates could not be formulated. The stability and error estimates de- 
veloped here for a method without an additional streamline diffusion term are made 
possible by: (1) the different and transparent format in which the problem is posed; 
and (2) the different machinery that is brought to bear on the problem. With respect 
to the last point, advantage is taken of some developments formulated by Ern and 
Guermond [7]. 

In the remainder of this work, the equation of interest and the numerical method 
to be analysed are first formalised. This is followed by analysis of the hyperbolic 
case, for which satisfaction of an inf-sup is demonstrated. The the diffusive limit 
case is then considered, for which demonstration of coercivity suffices. The results of 
some numerical simulations are then presented in support of the analysis, after which 
conclusions are drawn. 

2. Interface stabilised method. 

2.1. Model problem. Consider a polygonal domain ft C M rf , where 1 < d < 3, 
with boundary T = dQ. The unit outward normal vector to the domain is denoted 
by n. The advection-diffusion-reaction equation reads: 

(iu + a ■ Vu — kV 2 u = / in Q, (2-1) 

where [i > and K > are assumed to be constant, a : ft — > M. d is a divergence- 
free vector field that is Lipschitz continuous on fl and satisfies |o|| < 1, and 
/ : O — > M is a suitably regular source term. The divergence-free condition on a can 
easily be relaxed to /i — (1/2) V • a > 0. Portions of the boundary on which a ■ n > 
are denoted by T + , and portions on which a ■ n < are denoted by T_. A function 
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( is defined on boundaries such that ( = on outflow portions of the boundary (r + ) 
and C = 1 on inflow portions of the boundary (T_). 

For the case k > 0, the boundary is partitioned into IV and such that 
r^v U To = r and Tjy CiTd = 0, and the boundary conditions 

(—Cua + kVu) ■ n = q on Tn, 

y S ' B 2.2 

u = on r D , 

are considered, where g : IV — > R is a suitably smooth prescribed function. For the 
case k = 0, then To = 0, Tn = T_ and the considered boundary condition reads: 

— ua ■ n = g on r_. (2-3) 

2.2. The method. Let T be a triangulation of SI into non-overlapping simplices 
such that T = {K}. A simplex K e T will be referred to as a cell and a measure of 
the size of a cell K will be denoted by hx, with the usual assumption that Kk < 1, 
and /i = ma,XKeT hx- The boundary of a cell if is denoted by dK and the outward 
unit normal to a cell is denoted by n. The outflow portion of a cell boundary is 
the portion on which a ■ n > 0, and is denoted by dK + . The inflow portion of a 
cell boundary is the portion on which a n < 0, and is denoted by dK-. As for 
the exterior boundary, the function ( is defined such that ( — on dK + and £ = 1 
on dK-. The set of all facets T = {F} contained in the mesh will be used, as will 
the union of all facets, which is denoted by T°. Adjacent cells are considered to share 
a common facet F. 

The bilinear and linear forms for the advection-diffusion-reaction equation are 
now introduced. Using the notation w = (w, w) and v — (v, v), consider the bilinear 
form: 

B(w,v)= / fiwvdx+ / (— aw + kVw) ■ Vw dx 
Jn Jn 

+ / [ —aw + kVw — ( (a — - — n ) (w — w) | • n (v — v) ds 
^JokK V h K J K ') 



+ j k{w — w)Vv ■ nds + I a-nvwds (2.4) 
K JdK Jr + 

and the linear form 

L(v)= / fvdx+ / gvds, (2.5) 

where a > 0. The relevant finite element function spaces for the problem which will 
be considered read 

W h = {w h e L 2 (SI) , w h \ K e P k (K) VJfeT}, (2.6) 
Wh = {w h € H l (r°) ,w h \ F e P k {F)VF e T,w h = on T D } , (2.7) 

where < I < 1 and Pk (K) denotes the space of standard Lagrange polynomial func- 
tions of order k on cell K. The space Wh is the usual space commonly associated 
with discontinuous Galcrkin methods, and the space Wh contains Lagrange polyno- 
mial shape functions that 'live' only on cell facets and are single- valued on facets. The 
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choice of I, which determines the regularity of the facet functions at cell vertices in two 
dimensions and across cell edges in three dimensions, will have a significant impact 
on the structure of the resulting matrix problem. Using the notation Wt = Wh x Wh 
and Vh — (vh,Vh), the finite element problem of interest reads: find Uh € W£ such 
that 

B (u h , v h ) = L (v h ) Vv h e Wjf. (2.8) 

To motivate the terms appearing in the bilinear form, it is useful to consider 
the case in which v = and the case in which v = separately. Considering first 
Vh = (vh, 0), the variational problem corresponding to equation (|2.8[) for a single cell 
reads: given Uh € Wh, find ut & Wh such that for all Vh € Wh 

/ fJ,u h v h dx+ la (uh) ■ Vv h dx - / a (u h ) ■ nv h ds 

JK JK JdK 

n{u h - Uh)Vvh ■ nds = / fv h dx, (2.9) 

8K JK 

where a (w) — —a\7w + kVto is the usual flux vector and er (w) is a 'numerical flux', 

(a — - — n ) (w — w) . (2-10) 

h K J 

The problem in equation (|2.9[) is essentially a cell-wise postulation of a Galerkin 
problem for equation (I2.ip subject to the weak satisfaction of the boundary condition 
Uh = Uh- In the numerical flux, the presence of the term £ provides for upwinding of 
the advective part of the flux, and the term (oik/Iik) n{w — w) is an interior penalty- 
type contribution to the numerical flux The term J gK k [Uh — Uh) Vu/j - nds is 
typical of discontinuous Galerkin methods for elliptic problems, and resembles that 
in Arnold et al. @ for the Poisson equation. The numerical flux can be evaluated 
on both sides of a facet. On the outflow (upwind) portion of a cell boundary, the 
advective part of the numerical flux is equal to the regular advective flux. On the 
inflow (downwind) portion of a cell boundary, the advective part of the numerical flux 
depends on the interface function, taking on — au. The diffusive numerical flux on a 
cell boundary has contributions from the regular flux and a penalty-like contribution 
which depends on the difference between Wh and the interface function Wh- Setting 
Vh = 1 in equation (|2.9[) . 

/ fiuhdx— / &(uh)-nds= / f dx, (2-11) 

JK JdK JK 

which demonstrates local conservation in terms of the numerical flux. Note that the 
numerical flux defined in equation (|2.10[) is not single-valued on cell facets. Setting 
Vh = (0,v~h) furnishes the problem: given Uh S Wh, find ilh € Wh such that for all 

v h e w h 



E 

K 



I a (uh) ■ nvh ds + / a-nuhVhds~ / gvhds, (2-12) 
JdK Jr + Jr N 



which is a statement of weak continuity of the numerical flux across cell facets. 

Noteworthy in the bilinear form is that the functions Wh, which are discontinuous 
across cell facets, are not linked directly across facets. They are only linked implicitly 



AN INTERFACE STABILISED FINITE ELEMENT METHOD 



5 



through their interaction with uih- Setting Vh = (^,0) leads to a local (cell-wise) 
problem, which, given Uh and / can be solved locally to eliminate Uh in favour of Uh ■ 
This process is commonly referred to as static condensation. Then, setting Vh = 
(0,Vh), one can solve a global problem to yield the interface solution Uh- The field Uh 
can then be recovered trivially element- wise. To formulate a global problem with the 
same number of degrees as a continuous finite element method, I in equation (I2.7[) must 
be chosen such that there is only one degree of freedom at a given point; the interface 
functions are continuous at cell vertices in two dimensions and along cell edges in three 
dimensions. Further details on the formulation of the interface stabilised method and 
various algorithmic details can be found in Labeur and Wells Q. 

The formulation of Hughes et al. [1| can be manipulated into framework presented 
in this section, and in the hyperbolic limit coincides with the formulation presented 
here. In the case of diffusion, Hughes et al. [l| adopted an upwinded diffusive flux 
whereas the diffusive flux is centred in the present method. The formulation presented 
in Cockburn et al. 3] follows the same framework as Labeur and Wells Q, although 
the use of functions lying in 1? (r°) on facets is advocated. 

The method is now shown to be consistent with equation (12.11) . If u solves equa- 
tion (|2.ip . it is chosen to define u = (u,u). The action of the trace operator in the 
second slot is implicit in this definition (this will be expanded upon in Section [3]) . 
With this definition of u consistency can be addressed. 

Lemma 2.1 (consistency). If u = (u,u), where u € H m (f2) is a solution to (|2.1j) 
with m = 2 if k > and m — 1 otherwise, and if solves (|2.8|l . then for all v% € 

B(u-u h ,v h )=0. (2.13) 



Proof. Since Uh is a solution to (|2.8p and due to the bilinear nature of B, it 
suffices to demonstrate that B (u, Vh) — L (vh) = 0. Considering first B (u, (v^, 0)) — 
L ((vh,0)), which is presented in equation (I2.9[) . after applying integration by parts 

B (u, (v h , 0)) - L ({v h , 0)) = [ Uu + aVu - k\7 2 u - /) v h dx = 0, (2.14) 

since u satisfies (I2.1j) for k — 0. Considering now B (it, (0, Vh)) — L ((0, Vh))t which is 
presented in equation (|2.12l) . 

B(u,(0,v h ))-L((0,v h ))= [ ((-(ua + KVu)-n-g)v h ds = 0, (2.15) 

Jr N 

since u satisfies the boundary condition in (|2.3D . Summing equations (|2.14p and (|2.14j) 
and subtracting B (uh, Vh) — L (vh) = concludes the proof. □ 

2.3. Limit cases. The method will be analysed for the hyperbolic (k — 0) and 
elliptic (a = 0, \i = 0) limit cases. The bilinear form is therefore decomposed into 
advective and diffusive parts, 

B(w,v) = B A (w,v)+B D {w,v), (2.16) 

where 

Ba (w, v) = fiwv dx — / aw ■ \7v dx — 2. \ a ■ nw (^ — v ) ds 
Jn K Jk k Jok + 

— / a • nw (v — v) ds + / a ■ nwv ds (2.17) 
K J'dK_ Jr + 
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and 



Bd (w, v) = 2_] / kVw ■ Vv dx 

K Jk 

+ 'S~ y / ( kVw + - — n (w — w) ] ■ n(v — v) ds 

x JdK \ h K J 

+ V / k(w - w) Vv ■ nds. (2.18) 

K JdK 

Stability and error estimates will be proved by analysing Ba (w,v) and Bo (w,v) 
independently. 

2.4. Conventional discontinuous Galerkin methods as a special case. If 

the functions defined on facets are defined to be in L 2 fT ) (I = in equation (|2.7D ). 
then for the hyperbolic case the formulation reduces to the conventional discontinu- 
ous Galerkin formulation with full upwinding of the advective flux 

mini m. in the 



diffusive limit, it reduces to a method which closely resembles the symmetric interior 
penalty method II, HI ■ Of prime practical interest is the case where the interface 



functions are continuous as this leads to the fewest number of global degrees of free- 
dom, but the special case of I — is considered briefly in this section to illustrate a 
link with conventional discontinuous Galerkin methods. 

For the case ji = k = 0, setting Vh = everywhere and Vh = everywhere with 
the exception of one interior facet F, the method implies that at the facet F 

I aw h -n + v h ds= / aw h+ -n + v h ds, (2-19) 
Jf Jf_ 

where the subscript '+' indicates functions evaluated on the boundary of the upwind 
cell. This implies that for a given Wh, the facet function uih simply takes on the 
upwind value on each facet. Inserting this into equation (|2.17p and setting Vh = 0, 

B A (w h ,v h ) = / fiWhVh dx - J2 / aw h -Vv h dx 
Jn K Jk 

+ 2, \ a ' nw hVh ds + / o, ■ nwh+Vh ds, (2.20) 

K JdK + K JdK- 

which is the bilinear form associated with the classical discontinuous Galerkin formu- 
lation for hyperbolic problems with full upwinding. 

The diffusive case [k = 1, fx = 0, a = 0, a > 0) is now considered, in which case 
the subscripts '+' and '— ' indicate functions evaluated on opposite sides of a facet. 
Following the same process as for the hyperbolic case leads to 

-^-w h v h ds= - [ \-Vw h - ■ n_ + -^-Wh- ) Vh ds 
h K 2 J F \ h K J 

+ ~ J ^-Vw/ l+ • n + + J^p w h+^j Vh ds (2.21) 



on facets. Assuming for simplicity that hK is constant, inserting the expression for 
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Wh into (|2.18p and after some tedious manipulations, the bilinear forms reduces to: 



B D (w h ,v h ) / Vwft • \7v h dx - I (Vw h ) ■ lv h } ds 



(Vv h )ds+^-[ M • {v h J ds-^ [ lVw h JlVv h Jds, (2.22) 

where (a) = 1/2 (a + + a_) and [a] = (a + n + + a_n_) are the usual average and jump 
definitions, respectively. This bilinear form resembles closely that of the conventional 
symmetric interior penalty method, with the exception of the term which penalises 
jumps in the gradient of the solution. 

3. Notation and useful inequalities. The standard norm on the Sobolev 
space H S (K) will be denoted by ||-|| s K and the H S (K) semi-norm will be denoted by 
I 'L k- Constants c which are independent of Hk will be used extensively in the pre- 
sentation. The values of constants without subscripts may change at each appearance, 
and the value of any constant with a numeral subscript remains fixed. When c ap- 
pears with a parameter subscript, this indicates a dependence on a model parameter. 
For example, c M indicates a dependence on fi. 

Use will be made of various estimates for functions on finite clement cells for the 



case Iik < 1. In particular, use will be made of the trace inequalities [13l. Il4| 

II< 9 *<c(/%>|||U + M<k) VveH\K), (3.1) 
WVvnfcsK^c^MlK + hK \v\ 2 2tK ) VveH 2 (K). (3.2) 
On polynomial finite element spaces, the inverse estimate [l5, 14 1 
\ v h\i ,k 

<ch K || Vv h £P k (K) (3.3) 

will be used extensively. Combining equations (13.11) and (I3.3[) leads to 

Kllo.ax < eh'* \\v h \\ QtK Vv h GP k (K). (3.4) 

Frequently, functions defined on Q or on a finite element cell K will be restricted 
to an interior or exterior boundary. For finite element functions defined on a cell, 
owing to the continuity of the functions on a cell the trace is well-defined point-wise 
on the cell boundary. When considering functions in H s (fi) restricted to T°, the 
action of a trace operator 7 : H s (Q) — > H s ~ x / 2 (T°) should be taken as implied in 
the presentation. 

4. Analysis for the hyperbolic limit. The interface stabilised method is first 
analysed for the hyperbolic limit case which corresponds to the bilinear form in equa- 
tion (|2.17p . For this case the spaces 

W{h) = W h + H 1 (SI) , (4.1) 
W(h) = W h + H 1 ' 2 (r°) , (4.2) 

will be used in the analysis, as will the notation W*(h) = W(h) x W(h). The space 
W(h) has been defined such that it contains the trace of all functions in H 1 (SI) on T . 
This will prove important in developing error estimates. 
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Introducing the notation a n — \a-n\, two norms are defined on W*(h). The first 
is what will be referred to as the 'stability' norm, 



1 - 


i 


l 


2 


\\a n (v - v) 


+ 


a n v 




0,dK 




o,r 



IIMIIa = V IMIo.n + J2 h K \\a ■ Vv\\ 2 K + ^ 

K K 

The second norm, which will be referred to as the 'continuity' norm, reads 



\ v \\\a> = \\\ v \\\a + k K \\ V \\l,K + \\ a ™ V 



K 



K 



2 

0,Oif_ 



\\a n v 



K 



2 

0,dK + 



(4.3) 



(4.4) 



aftVh 



2 

O.dK 



Control of Vh & W£ in terms of the 1 1 1 - 1 1 1 ^ norm also implies control of Iik 

due to the following proposition. 

Proposition 4.1. There exists a constant c > such that for all K £ T and for 
all v h e W£ 



h K IKIIo,9K ^ C (iK ~ v h\\l. dK + hh\\l tK ) 



(4.5) 



Proof. Using the triangle inequality and the inverse inequality (|3.4I) : 

h K \\vh\\l idK =h K \\v h - v h + v h \\l >gK 

<h K (\\v h - V h \\ 0tdK + \\Vh\\o,dK) 



<2h K (\\v h -v h \\ 2 dK + chj/ \\v h \\ 2 Q K j 
<c(\\v h -v h \\l tdK + IKHo, 



(4.6) 



□ 

4.1. Stability. Stability of the interface stabilised method for hyperbolic prob- 
lems will be demonstrated through satisfaction of the inf-sup condition. Before consid- 
ering the inf-sup stability, a number of intermediate results are presented. The analysis 
borrows from the approach of Ern and Guermond Q to discontinuous Galerkin meth- 
ods (see also Ern and Guermond [14], Section 5.6]). A similar approach is adopted by 
Burman and Stamm [4]. 

Lemma 4.2 (coercivity) . For all v £ W*(h) 



B A (v,v) > HMIo,: 



2 ^ 

K 



1 - 


2 1 
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2 


\\a n (v - v) 




a n v 




OMR 2 




o,r 



(4.7) 



Proof. From the definition of P>a (v,v) and the fact that a is divergence-free, it 
follows from the application of integration by parts to (|2.17[) and some straightforward 
manipulations that 



1 

2 ^ 

K 



II - 


2 1 


l 


\\a n (v - v) 






0,dK 2 





2 

o,r 



(4.8) 



□ 
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As is usual for advection-reaction problems, B A (v, v) is coercive with respect to 
a particular norm, but the norm offers no control over derivatives of the solution. 
Consider a function Zh which depends on Wh £ W£ according to 

z h = (z h ,0) = (-h K a K -Vw h ,0), (4.9) 

where ax is the average of a on cell K. Lipschitz continuity of a implies the following 
bound on a cell K HIB|: 

\\a - a K \\ Loo(K) < ch K \a\ WL(K) . (4.10) 



Lemma 4.3. If the function Zh depends on Wh according to equation (14. 9[) . then 
for all Wh G Wt there exists a c\ > such that if Vh = c%Wh + Zh, then 



1 

2 Ill'WftllU < b a (z h ,w h ) + ciB A (w h ,w h ) = B A {v h ,w h ) . (4.11) 



Proof. Consider first two bounds on Using equation (|4.10D and the inverse 

estimate (|3.3[) . 



(4.12) 



W Z h\\o,K = \\ h Ka K ■ VWh\\ 0>K 

< \\h K a ■ Vw h \\ 0>K + \\h K (a - a K ) ■ Vw h \\ K 

< \\h K a ■ Vw h \\ K + ch K \a\ w ^ {K) \\h K Vw h \\ aK 

< \\h K a ■ Vw h \\ Q K + c\a\ WL(K) \\h K w h \\ 0K , 

and from the inverse estimate (|3.3[) 

\\ z h\\a.K = \\ h Ka K ■ S7w h \\ Q K < c\\a\\ Lao{K) \\w h \\ 0fK < c\\w h \\ K . (4.13) 
From the definition of the bilinear form in equation (|2.17l) , 



y^/tic ||a ■ Van ||o,k = B A (*h, ™h) + ^2 h K ^a K ■ (V 

K K 

+ / Jfeg / {a-Vwh)(a-a) K -Vwhdx 

is J K 



w h ) w h dx 



~y2 h K a n (w h - w h ) a K ■ Vwh ds. (4.14) 
K J 9K+ 

Applying the Cauchy-Schwarz inequality to the various terms on the right-hand side, 
h K \\a ■ Vw h \\l K < B A {z h , w h ) + ]T || 

K K 

+ ) \\hta ■ \7wh\\ \\hj c (a-a) K -'Vwh 



K 



0,K 



\\h K a K ■ ^w h \\ om+ \\a n (w h - w h )\\ o dK+ . (4.15) 



K 
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Each term is now appropriately bounded. Using equation (|4.13[) . 

llM^llo,x \\ h Ka K ■ Vw h \\ 0jK <c^Y^ \\ w h\\l lK ■ 



(4.16) 



K 



K 



Setting i?2 = Kix a ' ^ w h 



0,K 



tlx (a — ax) ■ VwJ and using (14.101) . an inverse 



\0,K 



inequality and Young's inequality, 



R 2 <c 



h^a ■ Vwh 



\ a - a K\\ L ^( K ) 



h 2 K Vw h 



0,K 



<ch K \ a \w^(K) h 2 K a-Vw h ^ h 2 K Vw h 



(2 



0,K 
2 



(4.17) 



<ch K \a\ WL(K) I — \\a ■ Vwh\\ 0>K + y ||w/i|| 0)J r 
where £2 > but is otherwise arbitrary. Setting £2 = 2c \a\ w i ^ 



(4.18) 



R2 < -^ik \\a ■ S7w h \\ Q K + ch K |a| w i_ (K) ||^|| 0)K . 

Setting i? 3 = \\h K a K ■ Vwh\\ 0>dK+ \\a„ (w h ~ Wh)\\ o gK+ and using equation (|4.12jl and 
Young's inequality, 

_ 1 

R3 <ch K 2 \\h K a K ■ Vw h \\ 0K \\a n {w h - w h )\\ o dK+ 



<c 



a 



<- 



< 



2e 3 
ch K 

£3 



0,K 



\W^(K) 



h 2 K w h 
h 2 K w h 



J \\a n {w h - w h )\\ 0idK+ 



O.K 



+ ^J-\Wn {Wh ~ W h )f Q<dK+ 



\\a ■ Vw h \\l tK + ( k) IWIo.k) + ~Y IK (™h - Wh)\\l tdK+ , 



(4.19) 



where £3 > but is otherwise arbitrary. Setting £3 = 4c, 



1 2 
R3 < -^ik \\a ■ Vwh\\ 0iK 



+ c\h K \a\ w ^ (K) \\w h \\ QiK 
Combining these results leads to 

\^ h K || a • ^w h \\l K < B A (z h ,w h ) 

ha l a lwi(n)\ 



lL°°(f2) 



an (w h - Wh) 



OMK, 



(4.20) 



K 



cmax la 



L°°(Q) 



Mlhhiio^ + Xllr" i w h-wh 

\ K 



0,8K 



(4.21) 

From the above result, the definition of the norm in (14.3[) and coercivity (|4.7|) . the 
lemma follows straightforwardly with ci = cmax (l, |a|^i ^ /fij. □ 

Proposition 4.4. For Zh which depends on Wh according to equation (|4.9I) . 
there exists a c-i > such that for all Wh £ Wt 

\\\zh\\\ A <c 2 \\\w h \\\ A . (4.22) 
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Proof. The components of HI^TilH^ can be bounded term-by-term. Using equa- 
tion gins}, 

A* \\zhfa.K = A* \\h K a K ■ S7w h \\ 2 QK < cfi \\a\\ 2 Lx{n) \\w h \\ 2 QK . (4.23) 
Using the inverse inequality f|3.3[) and equation (|4.12[) . 



h K \\a ■ Vzh\\ * < ch K \\ 



< c (h K \\a ■ Vw h \\l K + h K \a\ 2 w ^ [K) \\w h \\ 

For the facet term, 



2 

0,K 



(4.24) 



2 

anh K a K ■ Vw h 



0,dK 



< ch K \\a n \\ L oo {dK) (h K \\a-Vw h \\ 2 Q K + h K \a\ 2 w ^ (K) IWIq^J ■ (4.25) 

This proves that < cmax M|a|| ioo( - n ) , ^/i|a|^ (f2) / 'fij \\\w h \\\ A , with c 2 = 

cmax (ll«lli-(n) :l a l^(n) D 

Setting Vh = c\Wh + Zh, the preceding proposition also implies that 

IIKUL = HlcitOfc + z h \\\ A < (ci + c 2 ) \\\w h \\\ A . (4.26) 

Now, using the preceding two results, the demonstration of inf-sup stability is straight- 
forward. 

Lemma 4.5 (inf-sup stability). There exists a f3 A > 0, which is independent of 
h, such that for all Vh € Wj* 

B A (v h ,w h ) 

SU P -iii iii > PaWIvhWIa- (4.27) 

w h ew£ \\\Wh\\\ A 



Proof. For non-trivial — ciWh + Zh, combining Lemma I4.3l and Proposition 
(see also equation (|4.26p ) yields 

IIKIIU III^HL ^ ( c i + c 2) \\\ w h\\\ A < 2 (ci + c 2 ) B A (v h , w h ) , (4.28) 

which implies that for f3 A = 1/ (2 (ci + C2)), there exists a function G W 7 ^ such 
that 

^IKIIU^ ^y V« fc EWJf. (4.29) 

This is satisfaction of the inf-sup condition. □ 

Note the dependence of ft A on the problem data; it becomes smaller as gradients 
in a become large and as fi becomes small. In practice, this is a rather pessimistic 
scenario since often additional L 2 control will be provided by the prescription of the 
solution at inflow boundaries. Numerical experiments with /i = are usually observed 
to be stable. 
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4.2. Error analysis. To reach an error estimate, continuity of the bilinear form 
with respect to the norms defined in equations (|4.3p and (I4.4|) is required. It is 
the continuity requirement which necessitates the introduction of the norm |||-|||^, in 
addition to the stability norm HI'IH^. 

Lemma 4.6 (continuity). There exists a Ca > 0, which is independent ofh, such 
that for all w G W*(h) and for all Vh G 



\B A (w,v h )\<C A \\\w\\\ A , \\\v h \\\ A 
Proof. From the definition of the bilinear form: 



(4.30) 



\B A (w,v h 



/ pwj, dx — / aw ■ Vvh dx — \ . I a n w (vh — Vh) ds 
Jn K Jk k JdK + 



a n w (vh - Vh) ds + a n wv h ds 
OK- Jr + 



< 



+ E 

K 

Yl IH!o,K \\Vh\\ ,K + W a ■ V ^llo. 
K 



K 



K 



K 

II - 
\\anW 



1 1 




i 

an (Vh 


- Vh) 






0,dK+ 


0,dK + 


i - 




1 

a n {vh 


- v h ) 






0,dK- 


0,dK_ 



lo,r 4 



a n v h 



o,r 4 



Now, bounding each term, 



$^IHIo,K \\Vh\\ ,K 



■ < 



K 



EV \\ W \\o,K h K 11° • V v h\\()M ^ 
K 



K 



K 







an (Vh 


- v h ) 








II - 




1 








an {Vh 


- Vh) 




0,dK- 



< 



Q.dK, 



< 



0.8K- 



a n w 



o,r. 



aftVh 



< 



o,r 4 



1 w \ 


\a> 


Vh\ 


1 w \ 


\a> 


Vh\ 


1 w \ 


\a> 


Vh\ 


1 w \ 


\a> 


Vh\ 


1 w \ 


\a> 


Vh\ 



IA > 



(4.31) 

(4.32) 
(4.33) 
(4.34) 
(4.35) 
(4.36) 



Summation of these bounds leads to the result, and demonstrates that Ca = 1. □ 
The necessary results are now in place in to prove convergence of the method. 
Lemma 4.7 (convergence). For the case k — 0, if u — (u,u), where u solves 

equation (|2.1[) and Uh is the solution to the finite element problem (|2.8p . then 



\\u - u h \ 



< 1 



7T") ™i \\\u-v h \\\ A , 
Pa J v h ew* 



(4.37) 
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Proof. From inf-sup stability fLemma I4.5[) . consistency (Lemma 12.11) and conti- 
nuity of the bilinear form (Lemma I4.6[) : 

Q in in / B A (u h ~w h ,v h ) B A (u-w h ,v h ) 
PA\\\u h -w h \\\ A < sup - — = sup [ii [ii 

v h €W? \\\ v h\\\ A v h EW* \\\ v h\\\ A 

(4-38) 

^ n \\\u-w h \\ A , \\\v h \\ A 

<C A sup — f =CA\\\u-w h \\\ A , . 

v h ew* \\\ v h\\\ A 

Application of the triangle inequality 

\\\u - u h \\\ A < \\\u - w h \\\ A + \\\w h - u h \\\ A (4.39) 

and |||i>|IIa ^ IMIIa' yields the result. □ 

Lemma 4.8 (best approximation). For the case k — 0, if u G H k+1 (f2) solves 
equation (|2.ip and u = (u,u), and it/,, is the solution to the finite element prob- 
lem (|2.8p . then there exists a c^ a > such that 



and 



\u - u h \\\ A < c^. a h k+1 i \\u\\ k+l<n (4.40) 



\\u - u h || 0i n < c^ a h k+ ? \\u\\ k+i n . (4.41) 



Proof. The continuous interpolant of u is denoted by T^u = (l^u,I^u), where 
IhU E Wh n C (ft) and X^u — I/ju| r0 , which is contained in Wh- The standard 
interpolation estimate reads: 

\\u-l h u\\ nhK <ch h K +1 - m \u\ k+hK . (4.42) 

Bounding each term in ||it — I^ulH^, , 



\\u 


-Ihu\\l lK 


< ch 2{k+1) \u\ 2 

— UL K \ u \k+l,K ) 


(4.43) 


hx ||a ■ V (u - 


-^ u )llo,x 


ul K \ a \k+l,K J 


(4.44) 


u — IhUj — (u — 






(4.45) 


h K X \\u 




< rh 2k+1 \v\ 2 

ul K \ a \k+l,K ) 


(4.46) 



\\u -l h u\\ QdK = \\u- 1 h u\\ odK 

< c ( h K \\u -l h u\\ 0K + h K \u -l h u\ l K \ < ch K + \u\ k+1K . (4.47) 

Using these results and equation (|4.37l) leads to the convergence estimates. □ 

5. Analysis in the diffusive limit. The diffusive limit (a = 0, fj, = 0) is now 

considered, in which case the bilinear form is given by equation (|2.18l) . The analysis 
of the diffusive case is considerably simpler than for the hyperbolic case since stability 
can be demonstrated via coercivity of the bilinear form. Analysis tools and results 
which are typically used in the analysis of discontinuous Galerkin methods for elliptic 
problems [8[ are leveraged against this problem. 
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To ease the notational burden, the case of homogeneous Dirichlet boundary con- 
ditions on r is considered. The extended function spaces 



W(h) = W h + H 2 (Cl) n Hi (ct) , 

w{h) = w h + h 3 q /2 (r°) , 



(5.1) 
(5.2) 



will be used, where H^ 2 (T°) denotes the trace space of H 2 (Cl)nHQ (Cl) on facets L°. 

As for the hyperbolic case, two norms on W* (h) = W(h) x W(h) are introduced 
for the examination of stability and continuity. The 'stability' norm reads 



K 

and the 'continuity' norm reads 



K 



UK 
h~ 



K 



\o,dK ' 



\D 



E 

K 



hxK . ,2 

a 1 h - K 



(5.3) 



(5.4) 



It is clear from the definitions that |||i>|||^) < HMH^,, , but there also exists a constant 
c > such that for all Vh € W£ 



UKIIU<c(i + , 

since from equation (|3 . 1 [) it follows that 



D • 



h 2 K \ v h\ 2 K ^ c ll V ^lloF/ VVhEWh- 



(5.5) 



(5.6) 



Therefore, the norms 

To demonstrate that 
facet F 



D 



and 



D 



and 



y are equivalent on the finite element space Wt. 
, do constitute norms, first recall that for a 



D 



v -\\o,F =ElK W + 
F 

<£n* 



I0,F 



+ II0.F 



+ U — V- 



0,F 



(5.7) 



'Eli- 



V 



O.dK 



K 



1 




K 2 V 


< Cl ( 







Denoting the average size of two cells sharing a facet by /if, 

on - Cl ( E K ll Vy llo,K + E IT" H U+ " Ho,F ) " Cl t 1 + ' 

V K F F / 

(5.8) 

where the first inequality is a standard result (See Arnold |13l . Lemma 2.1] and Ern 
and Guermond 14|, Lemma 3.45]). Hence, \\\-\\\ D and constitute norms. 

5.1. Stability. Before proceeding to coercivity of the bilinear form, an interme- 
diate result is presented. 

Proposition 5.1. There exists a constant c > such that for any e > and all 

v h e w* 



IdK 



kWvh -n(v h -v h ) ds 



< ecu \\Vv h f 0tK + -t— |K - v h \\l. 



dK ■ 



(5.9) 
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Proof. Applying to the term Jg K ^Vf,. ■ n(vh — Uft) ds the Cauchy-Schwarz in- 
equality, the inverse estimates and Young's inequality, 



OK 



nX7v h ■ (v h -v h )nds 



<2nh K \\Vvh ■ n\\ ogK \\v h - v h \\ Q 



OK 



<h K K€ \\Vv h ■ nf Q dK + — \\v h - v h \\l ok 
eriK 



<ecn 



(KIi.jc 



+ h 2 K \v h \l K ) +-7— IK -Ufello, 
/ en K 

K II- l|2 



OK 



<ecK \\Vv h \\ K + — \\v h - v h \\ Q dK , 



(5.10) 



which complete the proof. □ 

Lemma 5.2 (coercivity) . There exists a (3d > 0, independent ofh, and a constant 
ao > such that for a > ao and for all Vh 6 Wfc 



B D (v h ,v h ) > 13 D \\\v h \\\ D 

and there exists an ct\ > a?o such that for all S W£ 

1 2 
B D {vh,v h ) > - \\\v h \\\ D . 



(5.11) 



(5.12) 



Proof. Setting Wh = Vh in the bilinear form for the diffusive limit case (|2.18|) 
B D (v, v)=y^n \\Vv h |1 o K + 2 V] / K ^ v h -n(v h -v h ) ds 

J OK 



OK 

t— \\vh ~ v h \\ o dK . (5.13) 



h 

K 

Using Proposition 15 . 1 1 to bound the term J2k Iok K ^ v h ' n {^h — Vh) ds, 

B D (v h ,v h )>J2( 1 -^)K\\Wv h \\l K + ^(a--)^\\v h -v h \\l dK . (5.14) 



K K v ' 

Setting e = 1/Sc, where 6 > 1 but is otherwise arbitrary, 



K 



B D {v h ,v h ) > ( 1 - T ) «l|V^||o )if + ~ 5c) \\v h -v h \\l dK . (5.15) 



k x if " 

This proves equation (|5.11[) and demonstrates that «o = c. Setting (5 = 2, 

(Vh,V h ) > Y 2 K H V ^Ho,/f + H (« - 2c ) J— IK - V h\\o,dK ' ( 5 ' 16 ) 
K K K 

which proves (|5 . 1 2[) when a\ = 1/2 + 2c. □ 

The proof to Lemma 15.21 demonstrates that stability is enhanced for a larger 
penalty parameter. Stability demands that a > c, and when this is satisfied 8 can be 
chosen such that (3d approaches zero as a approaches ao, and such that /3d approaches 
one as a becomes much larger than ao- 
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5.2. Error analysis. The error analysis proceeds in a straightforward manner 
now that the stability result is in place. 

Lemma 5.3 (continuity). There exists a Co > 0, independent of h, such that for 
all w G W*(h) and for all v E W*(h) 

\B D (w,v)\<C D \\\w\\\ D ,\\\v\\\ D ,. (5.17) 

Proof. From the definition of the bilinear form, 
\B D (w, v h ) | < K \\Vw\\ Q K \\Vv h \\ Q K + II Vw ■ n\\ o dK \\v h - v h \\ Q dK 

K K 

+ ^2k\\w- w\\ o dK \\Vv h ■ n\\ o dK + 7— \\w - w\\ 0ydK \\v h - v h \\ 0t Q K ■ ( 5 - 18 ) 

K K K 

Each term can be bounded appropriately, 

J2 K W Vw Wo,K l|V»h|| 0iX < |||t»||| D , J2 ^ W^Wo,K > (5-19) 

K K 

^k\\Vw ■ n\\ QdK \\v h - v h \\ odK <^c(/i K 5 \w\i,k + h K \ w \2,k) \\vh ~ v h \\ 0t9K 

K K 

l 

<c^max(l,a^) |H|| B , £ ( \\v h - v h \\ om , (5.20) 

K K ^ nKj 

^Tk\\w- w\\ Q dK \\Vv h ■ n\\ Q dK < \\w - w\\ o dK (h K * \v\ 1K + h% |u| 2 ,jr) 

K K 

< cmax (l, aT*) \\\w\\\ D , £ (^) 2 (\v\ 1K + h K \v\ 2K ) , (5.21) 



J2~h^\\ W ~ W hdK\\ V - V \\o,dK < \\\ w \\\d>J2 (j^J \\v- V \\ ,dK- ( 5 - 22 ) 

Summing these inequalities shows that the bilinear form is continuous with respect 
to IH'llljj/, with Ca = cmax ML, a~^j . □ 

The penalty term a is usually taken to be greater than one, in which case Cd = c. 

Lemma 5.4 (convergence). If a is chosen suitably large such that the bilinear 
form is coercive, then if u solves equation (|2.1D and u — (u, u), and Uh is the solution 
to equation (|2.8I) for the case fi = 0, a = and k > 0, then 

III*" - Uh\\\ D < \ 1 + o ^ CD ) in f III" - w h\\\ D , ■ (5.23) 

y Pd J whew* 

Proof. Using coercivity, consistency and continuity: 

@d \\\u h - w h \\\ D < B D (u h - w h ,u h - w h ) 

= B D (u - w h ,u h - w h ) (5-24) 
< C D \\\u - w h \\\ D , \\\u h - w h \\\ D , , 
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and then exploiting IH^IH^,, < (l + ca *) \\\vh\\\ D (see equation (|5.5|l ). 

UK - w h \\\ D < Pn l C D (1 + ca- 1 ) \\\u - w h \\\ D , , (5.25) 

which followed by the application of the triangle inequality yields the desired result. 
□ 

Lemma 5.5 (best approximation). For the case \i = 0, a = and k > 0, if 
u G H k+1 (f2) solves equation (|2.1[) and u = (u,u), and is the solution to the finite 
element problem (|2.8[) , and a is chosen such that the bilinear form is coercive, then 
there a exists a c a > such that 

\\\u - u h \\\ D < c a h k \\u\\ k+hn (5.26) 

and 

ll"-^l!o,Q< c ^ fc+1 IMI fc+ i^- ( 5 - 27 ) 

Proof. The first estimate follows directly from the standard interpolation estimate 
for the continuous interpolant T^u = (lhU,Ihuj , where again IhU € Wh H C (f2) and 
I/jtt = Ihu\ r o, which is an element of Wh- Applying the standard interpolation 
estimate (|4.42p to |||ti — Z^uHl^,, 

\\V(u-I h u)\\ 2 0K <ch 2k \u\l +1K , (5.28) 
|| (it — Thu) - (u -2 h u)\\ 2 Q dK = 0, (5.29) 
h% \u-lu h \ 2 2 K < ch 2k \u\ 2 k+1K . (5.30) 

Using these inequalities leads to equation (|5 . 26|) . The L 2 estimate follows from the 
usual duality arguments. Owing to adjoint consistency of the method (since the 
bilinear form is symmetric), if w S H 2 (il) DHq (ft) is the solution to the dual problem 

B D (v,w)= [ (u-u h )-vdx yv£W*(h), (5.31) 
Jn 

and wi G W£ is a suitable interpolant of w, then from consistency and continuity of 
the bilinear form, and the estimate in (|5 . 26[) . it follows that 

\\ u ~ u h\\l.n = B D (u- u h ,w) 

= B D (u-u h ,w -wi) ^ 

< C D \\\u - u h \\\ D , \\\w - wi\\ D , 

< c a h\\w\\ 2n \\\u - u h \\\ D , . 

Finally, using the elliptic regularity estimate ||i/;|| 2 < c a \\u — Uh\\ q leads to 

ll u ~ u h\\o,n ^ C a h W\ u ~ U h\\\ D > ■ (5.33) 
The L 2 error estimate follows trivially. □ 



18 



GARTH N. WELLS 



6. Observed stability and convergence properties. Some numerical exam- 
ples are now presented to examine stability and convergence properties of the method. 
In all examples, the interface functions are chosen to be continuous everywhere (I — 1 
in equation (|2.7[) ). so the number of global degrees of freedom is the same as for a 
continuous Galerkin method on the same mesh. When computing the error for cases 
using polynomial basis order fc, the source term and the exact solution are interpo- 
lated on the same mesh but using Lagrange elements of order k + 6. Likewise, if 
the field a does not come from a finite element space it is interpolated using order 
k + 6 Lagrange elements. Exact integration is performed for all terms. All meshes 
are uniform and the measure of the cell size hx is set to two times the circumradius 
of cell K. 

The computer code used for all examples in this section is freely available in 



the supporting material [17( under a GNU Public License. The necessary low-level 
computer code specific to this problem has been generated automatically from a high- 
level scripted input language using freely available tools from the FEniCS Project [18l 



|20L [2l|. The computer input resembles closely the mathematical notation and 

ntag e 

of automation developments for methods that involve facet integration 19]. 



abstractions used in this work to describe the method. Particular advantage is taken 

m 

6.1. Hyperbolic problem. Consider the domain ft = (—1,1), with /x 
a = (0.8, 0.6), k = and u = 1 on T_. The source term / is chosen such that 



u = 1 



+ sin(7r(l + x)(l + y) 2 /8) (6.1) 



is the analytical solution to equation (|2.ip . This example has been considered previ- 
ously for discontinuous Galerkin methods by Bey and Oden [22[ and Houston et al. 
& 

The computed error \\\u — it/ilH^ is presented in Figure RPl for /i-refinement with 
various polynomial orders. As predicted by the analysis, the observed converge rate 
is k + 1/2. For all polynomial orders the method converges robustly. 

6.2. Elliptic problem. A problem on the domain = ( — 1, l) 2 is now consid- 
ered, with fi = 0, a = (0, 0) and k — 1. The source term / is selected such that 

u = sin (ttx) sin (ny) (6-2) 



is the analytical solution to equation (|2.1[) . The value of the penalty parameter is 
stated for each considered case. 

The computed errors in the L 2 norm for /i-refinement with elements of varying 
polynomial order and a — 5 are shown in Figure 16.21 In all cases, the predicted 
k + 1 order of convergence is observed. The computed results for a = 6 are shown in 
Figure 16.31 in which the convergence for the k — 2 case is somewhat erratic. Using 
a = 4fc 2 , since the penalty parameter for the interior penalty method usually needs 
to be increased with increasing polynomial order, reliable convergence behaviour at 
the predicted rate is recovered, as can be seen in Figure 

6.3. Advection-diffusion problems. An advection-diffusion problem is con- 
sidered on the domain ft = (— 1, 1) , with fi = 0, a = (e x (y cosy + sin y), e x y sin y) 
and for various values of n. The source term / is chosen such that equation (|6.2p is 
the analytical solution. For all cases, a = 4k 2 . 

The convergence behaviour is examined in terms of \\\u — ~\~ \\\u — ^/illl^)- 

The computed error for the case k = 1x 10~ 3 is presented in Figure 16.51 For this 
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advection dominated problem, the method is observed to converge at the rate k + 
1/2. For k = 0.1, the observed convergence response is presented in Figure WM A 
convergence rate of k is observed for the lower order polynomial cases, and the rate 
appears approach k + 1/2 for the higher-order polynomial cases. For k = 10, which is 
diffusion dominated, the observed convergence is presented in Figure [fTTl As expected, 
a convergence rate of k is observed for the diffusion-dominated case. 
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= 1 


H ^ ft 


= 2 




= 3 


' * ft 


= 4 




= 5 



10° 



Fig. 6.3. Convergence for the elliptic case in L 2 with h-refinement for various polynomial 
orders and a = 6. 






= 1 


H ^ jfe 


= 2 




= 3 


' * ft 


= 4 




= 5 



10" 



Fig. 6.4. Convergence for the elliptic case in I? with h-refinement for various polynomial 
orders and a = 4k 2 . 



7. Conclusions. Stability and error estimates have been developed for an in- 
terface stabilised finite element method that inherits features of both continuous and 
discontinuous Galcrkin methods. The analysis is for the hyperbolic and elliptic limit 
cases of the advection-diffusion-reaction equation. While the number of global de- 
grees of freedom on a given mesh for the method is the same as for a continuous finite 
element method, the stabilisation mechanism is the same as that present in upwinded 
discontinuous Galerkin methods. This is borne out in the stability analysis, which 
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Fig. 6.6. Convergence for the advection- diffusion problem in |||u — it/JII^ + W u ~ u hMo f or 
K = 0.1 with h-refinement for various combinations of k and a = 4k 2 . 



demands consideration of an inf-sup condition. Analysis of the method shows that it 
inherits the stability properties of discontinuous Galerkin methods, and that it con- 
verges in L 2 at a rate of k + 1/2 in the advective limit and k + 1 in the diffusive limit, 
as is typical for discontinuous Galerkin and appropriately constructed stabilised finite 
element methods. The analysis presented in this work provides a firm theoretical basis 
for the method to support the performance observed in simulations in other works. 
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Fig. 6.7. Convergence for the advection- diffusion problem in \\\u — Uh\\\ A + |||u — u/JH^ for 
k — 10 with h-refinement for various combinations of k and a = 4k 2 . 



The analysis results are supported by numerical examples which considered a range 
of polynomial order elements. 
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